#include <iostream>
#include <array>
#include "Core/TensorExpr.h"
#include "Embedded/EmbeddedGrid.h"
#include "Embedded/EmFuncFiller.h"
#include "Embedded/EmLevelOp.h"
#include "Embedded/Convection.h"
#include "../example/common/Ellipsoid.h"
#include "../example/common/TestUtility.h"

//const Real maxVel = 1.2;
const Real gTol = 1e-12;

template <int Dim, int desiredOrder, class T1, class T2>
std::array<Real,3> testCvt(int N, const T1 &exactVel, const T2 &exactCvt)
{
  // intialize the computational domain
  using rVec = Vec<Real, Dim>;
  const int numGhost = desiredOrder / 2;
  RectDomain<Dim> rd(Box<Dim>(0, N-1), 1.0/N, CellCentered, numGhost);
  Ellipsoid<Dim> ellipsoid({0.5, 0.5, 0.5}, {0.25, 0.125, 0.25});
  EmbeddedGrid<Dim> emGrid(rd, [&](rVec q, rVec &closest, rVec &normal, bool &contained, Real gTol) {
    ellipsoid.check(q, closest, normal, contained, gTol);
    normal = -normal;
    contained = !contained;
  }, true, 0.1, gTol);
#ifndef NDEBUG
  if(get_dbglevel() >= 2) {
    emGrid.dump(tmpos);
    tmpos.flush();
  }
#endif
  // fill the initial values
  EmFuncFiller<Dim> eff(emGrid);
  EmLevelOp<Dim>    elv(emGrid);
  Box<Dim> gbx = rd.getGhostedBox();
  Tensor<Real, Dim+1> phis(gbx, Dim);
  Tensor<Real, Dim> phi[] { {phis, 0},
                            {phis, 1},
                            {phis, 2} };
  Tensor<Real, Dim-1> bcData(emGrid.bdryNodes.size());
  for(int d = 0; d < Dim; ++d) {
    auto vel_d = [&](const rVec &x) { return exactVel(x)[d]; };
    eff.fillPointValues(phi[d], CellCentered, vel_d);
//    eff.fillBdryValues(phi[d], "curvilinear", vel_d);
    eff.fillBdryValues(bcData, "curvilinear", vel_d);
    eff.transferBdryValues(phi[d], bcData);
//    elv.zeroExtNodes(phi[d], false);
  }
  Tensor<Real, Dim+1> errs(gbx, Dim);
  Tensor<Real, Dim> err[] { {errs, 0},
                            {errs, 1},
                            {errs, 2} };
  for(int d = 0; d < Dim; ++d)
    eff.fillPointValues(err[d], CellCentered, [&](const rVec &x) { return exactCvt(x)[d]; });
  // carry out the computation
  Convection<Dim, desiredOrder+1> cvtOp(emGrid, phis);
  Tensor<Real, Dim+1> ress(gbx, Dim);
  Tensor<Real, Dim> res[] { {ress, 0},
                            {ress, 1},
                            {ress, 2} };
  cvtOp.apply(phis, ress);
  elv.zeroExtNodes(ress);
  errs = ress - errs;
  elv.zeroExtNodes(errs);
#ifndef NDEBUG
  if(get_dbglevel() >= 2)
    errs.dump(tmpos);
#endif

  std::array<Real, 3> errnorm;
  for(int p = 0; p <= 2; ++p) {
    rVec v;
    for(int d = 0; d < Dim; ++d)
      v[d] = elv.NRM(err[d], CellCentered, p);
    errnorm[p] = norm(v, p);
  }
  return errnorm;
}

template <int Dim, int desiredOrder, class T1, class T2>
void doTest(const std::vector<int> &gridSize, const T1 &exactVel, const T2 &exactCvt)
{
  const int numGrid = gridSize.size();
  std::vector<std::array<Real, 3>> errnorm;
  for(int n=0; n<numGrid; ++n) {
    std::cout << "\n\nTesting 1/h = " << gridSize[n] << std::endl;
    errnorm.push_back(testCvt<Dim, desiredOrder>(gridSize[n], exactVel, exactCvt));
  }
  printConvergenceTable(&gridSize[0], errnorm);
}

int main()
{
  reset_dbglevel(1);
  std::cout << "SpaceDim = " << SpaceDim << std::endl;
  std::cout << "===================================================================================" << std::endl;
#if DIM == 2
  auto exactVel = [](const Vec<Real, 2> &x) {
    return Vec<Real, 2> { ipow<2>(sin(M_PI*x[0])) * sin(2*M_PI*x[1]),
                          -sin(2*M_PI*x[0]) * ipow<2>(sin(M_PI*x[1])) };
  };
  auto exactCvt = [](const Vec<Real, 2> &x) {
    return Vec<Real, 2> { ipow<2>(sin(M_PI*x[0])) * sin(2*M_PI*x[0]) * ipow<2>(sin(M_PI*x[1])),
                          ipow<2>(sin(M_PI*x[1])) * sin(2*M_PI*x[1]) * ipow<2>(sin(M_PI*x[0])) } * (2.0*M_PI);
  };
  std::cout << "Testing fourth-order ..." << std::endl;
  doTest<2, 4>({256, 128, 64, 32}, exactVel, exactCvt);
  std::cout << "Testing sixth-order ..." << std::endl;
  doTest<2, 6>({256, 128, 64, 32}, exactVel, exactCvt);
  std::cout << "Testing eighth-order ..." << std::endl;
  doTest<2, 8>({128, 64, 32}, exactVel, exactCvt);
#elif DIM == 3
  auto exactVel = [](const Vec<Real, 3> &x) {
      return Vec<Real, 3> { ipow<2>(sin(M_PI*x[0])) * sin(2*M_PI*x[1]) * sin(2*M_PI*x[2]) * (1.0/2),
                            sin(2*M_PI*x[0]) * ipow<2>(sin(M_PI*x[1])) * sin(2*M_PI*x[2]) * (1.0/2),
                            sin(2*M_PI*x[0]) * sin(2*M_PI*x[1]) * ipow<2>(sin(M_PI*x[2])) * (-1.0) };
  };
  auto exactCvt = [](const Vec<Real, 3> &pos) {
      const Real x = pos[0], y = pos[1], z = pos[2];
      Real r0 = (M_PI/2) * ipow<2>(sin(M_PI*x)) * (cos(M_PI*x) * sin(M_PI*x) * ipow<2>(sin(2*M_PI*y)) * ipow<2>(sin(2*M_PI*z))
          + cos(2*M_PI*y) * sin(2*M_PI*x) * ipow<2>(sin(M_PI*y)) * ipow<2>(sin(2*M_PI*z))
          - 2 * cos(2*M_PI*z) * sin(2*M_PI*x) * ipow<2>(sin(2*M_PI*y)) * ipow<2>(sin(M_PI*z)));
      Real r1 = (M_PI/2) * ipow<2>(sin(M_PI*y)) * (cos(2*M_PI*x) * ipow<2>(sin(M_PI*x)) * sin(2*M_PI*y) * ipow<2>(sin(2*M_PI*z))
          + cos(M_PI*y) * ipow<2>(sin(2*M_PI*x)) * sin(M_PI*y) * ipow<2>(sin(2*M_PI*z))
          - 2 * cos(2*M_PI*z) * ipow<2>(sin(2*M_PI*x)) * sin(2*M_PI*y) * ipow<2>(sin(M_PI*z)));
      Real r2 = (-M_PI) * ipow<2>(sin(M_PI*z)) * (cos(2*M_PI*x) * ipow<2>(sin(M_PI*x)) * ipow<2>(sin(2*M_PI*y)) * sin(2*M_PI*z)
          + cos(2*M_PI*y) * ipow<2>(sin(2*M_PI*x)) * ipow<2>(sin(M_PI*y)) * sin(2*M_PI*z)
          - 2 * cos(M_PI*z) * ipow<2>(sin(2*M_PI*x)) * ipow<2>(sin(2*M_PI*y)) * sin(M_PI*z));
      return Vec<Real, 3> {r0, r1, r2};
  };
  std::cout << "Testing fourth-order ..." << std::endl;
  doTest<3, 4>({128, 64, 32}, exactVel, exactCvt);
  std::cout << "Testing sixth-order ..." << std::endl;
  doTest<3, 6>({128, 64, 32}, exactVel, exactCvt);
  //  std::cout << "Testing eighth-order ..." << std::endl;
  //  doTest<3, 8>({128, 64, 32}, exactVel, exactCvt);
#endif
  return 0;
}
